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Abstract : This paper contains some definitions, abbreviations, acronyms, concepts and coding of 
Numerical Computation involving q method. Numerical Computation plays an imperative role in solving 
real time and real life problems of engineering, mathematics and physics. It is an approach for solving 
complex mathematical problems using arithmetic operations. This approach involves formulation of 
mathematical models of physical situations that can be solved with mathematical operations. Applications 
of computer oriented numerical methods have become an integral part of life of all the modern scientists. 
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I. Introduction 

The field of numerical analysis predates the discovery of modern computers by many centuries. Linear 
interpolation was already in use for more than ten centuries. Many great mathematicians of the past were 
engrossed in study of numerical analysis and it is also apparent from the names of important algorithms 
like Newton’s method, Lagrange interpolation polynomial, Gaussian elimination, or Euler’s method.To facilitate 
computations by hand, large books were produced with formulas and tables of data such as interpolation points 
and function coefficients. Using these tables, often calculated out to sixteen decimal places or more for some 
functions, one could look up values to plug into the formulas given and achieve very good numerical estimates 
of some functions. 


The mechanical calculator was also developed as anapparatus for hand computation. These calculators evolved 
into electronic computers in the first generation of computers, and it was then found that these computers were 
also useful for administrative purposes. The discovery of the computer also influenced the field of numerical 
analysis, since now longer and more complicated calculations could be done. Numerical computing with the 
help of some special functions enhanced computational techniques. 

Heine’s basic hypergeometric series, or hypergeo metric q- series, are ^-analogue generalizations 

of generalized hypergeometric series, and are in turn generalized by elliptic hypergeometric series. A series x n is 
called hypergeometric if the ratio of successive terms x n+ i/x n is a rational function of n. If the ratio of successive 
terms is a rational function of q n , then the series is called a basic hypergeometric series. The number q is called 
the base. The basic hypergeometric series was first considered by Eduard Heine (1846). It becomes the 
hypergeometric series F( a,p;y;x) in the limit when the base q is l.Some of the q analogues are explained below. 


1.1 q Analogue of Exponential Function 

#-exponentialis dibasic analogue [Exton(1983)] of the exponential function, namely the eigen function of a q- 
derivative. There are many g- derivatives, for example, the classical ^-derivative, the Askey-Wilson [Wilson et 
al. (1985)] operator, etc. Therefore, unlike the classical exponentials, ^-exponentials are not unique. Three 
variants of exponential functions are given below. Third one is generalized formula and we can get first and 
second by changing values of cs. 
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1.2 q Analogue of /toiclntegration : The inverse operation to basic differentiation has also been discussed at 
some length by F.H. Jackson [Jackson (1904)]. This is represented by the symbol f 0(x) d ( qx ) and is referred 

to as ^-integration or basic integration. When q tends to unity, the basic integral reduces to the ordinary integral. 
The operations of basic differentiation and integration correspond exactly in every way to ordinary 
differentiation and integration of which these are generalizations. 


= Cl - q){&£r=ofl 1 7(<J r iO - oEr=o<? r /(<J r o)} (1.4) 
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fix)d{q,x) = (1 - q)Sf=_ K q7(q L ) (1.7) 


1.3 q Analogue of Trigonometric Functions :As like exponential function, various analogues of circular 
functions have been introduced. Jackson gave a class of q- circular function [Exton (1983)], [Gasper et al. 
(2004)] .He also introduced a class of circular function from point of view of pseudo-periodicity. 
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1.4 Properties of Trigonometric Functions 

sin q (x) +- cos q (x)cos t ^(x) = 1 

(1-10) 


ros^Cx) cos 1 i ;.-(.x) — sin^ .v (x) = cos,. (2x3 
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1.5 Basic Differentiation operator: Jackson introduced the operative symbol [Exton (1983)], [Gasper et al. 
(2004)] for basic differentiation defined by the relation A{0 (x)} = {(x) - (qx)} x' 1 (l-qf 1 . The operation of basic 
differentiation is defined by [Jackson (1904)] the relations 


B qx 0(x) = 
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, where x and q may be real or complex. This becomes the same as ordinary differentiation as the base q tends to 
unity. In order to avoid the possibility of perplexity with the ordinary difference operator, we shall write 

instead of A. Furthermore, the subscripts q and x will be omitted provided that there is no chance of vagueness. 
It will be seen that the possibility now arises of the existence of certain types of difference equations based upon 
this operator. 

(li3) 


1.6 Basic analogue of Taylor’s Theorem Jackson introduced [Exton (1983)] q analogue of Taylor’s Theorem 
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1.7 Variants of Laplace Transform :Hahn [Hahn (1949)] defined two analogues of Laplace Transform by the 
help of the integral equations 
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1.8 Basic Analogue of Heine’s Series : This series can be represented as [Heine (1847)], [Exton (1983)], 
[Gasper et al. (2004)] 
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1.9q Analogue of Euler’s Identity 
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1.10 Heine Equation :The Gauss hyper geometric function [Exton (1983)], [Gasper et al. (2004)] 2F1 (a, b; c; 
x) is a particular solution of the equation 

x(l — x)y" f [c - (l f c I b)x}y f E aby = 0 (1.21 a) 

which may be written in operational form as 

x(S E t0(6 -f iOy - 6(6' + c - l)y = 0. (1.21 b) 

If we replace the symbolic operations by their basic analogues, we obtain the ^-differential equation 
xl6 -f c: qilS + i>.; £g]y - [6; g][6 + c- 1; g]y = 0 (1.22) 

which on expansion takes the form 

x{q € - q a+k+1 x}B 2 y + He; g] - (g c [l E fa: g] E g & U.; g]M - [o; g][fc g]y = 0 (1.23) 


This is one of an infinite number of possible g-analogues of the hyper-geometric equation. 
l.llg-Gauss [Gasper et al. (2004)] summation formula :Gauss summation formula can be described as 
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1.12g-Plaff-Saalschutz’s summation [Gasper et al. (2004)] formula 

It can be described by the formula 


(1.24) 
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1.13 Some identities of g-shifted factorials [Exton (1983), Gasper et al. (2004)] are 
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1.14</-theta function 

g-theta function is given by 

0k;q) = n= =0 Cl (1.29) 

It can also be expressed as 

&0c\q) = Or; q)* (1.30) 


1.15 Hahn-Exton ^-Bessel function :Hahn-Exton g-Bessel function or the third Jackson g-Bessel function is 
a ^-analogue of the Bessel function, introduced by Hahn [Hahn et al. (1953)] in a special case and 
by Exton [Exton (1983)] in general. 


The Hahn-Exton g-Bessel function is given by 
v ' ^ HaEJ 


(1-31) 


1.16</-Weibull distribution :It is a probability distribution that generalizes the Weibull distribution [Gasper et 
al. (2004)] and the Lomax distribution (Pareto Type II). It is one example of a Tsallis distribution. 
The probability density function of a g- Weibull random variable is 


/U; q,A, ftj = 


0 x < 0 


, (—f 


, x>=0; 


(1.32) 


where g < 2, k> 0 are shape parameters and X > 0 is the scale parameter of the distribution. 

1.17Jackson ^-Bessel function (or basic Bessel function) :It is one of the three ^-analogues of the Bessel 
function introduced byF.H. Jackson [Jackson (1905)]. The third Jackson g-Bessel function is the same as 
the Hahn-Exton g-Bessel function. 

The three Jackson g-Bessel functions [Exton (1983)] are given in terms of the Pochhammer symbol and 
the basic hypergeometric function cp by 
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(1.34) 

(1.35) 


1.18Hopf algebra :Hopf algebra, named after Heinz Hopf [Hopf (1941)]. The representation theory of this 
algebra is predominantly fine, since the existence of compatible co -multiplication, co-unit, and antipode allows 
for the construction of tensor products of representations, minor representations as well as dual representations. 

1.19Quantum Affine Algebra (Affine Quantum Group) :It is a Hopf algebra [Drinfeld (1985)] that is a q- 
deformation of the enveloping algebra of an affine Lie algebra. It was introduced by Drinfeld and Jimbo. It a 
particular case of their normal construction of a quantum group from a Cartan matrix. 
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1.20 Hecke Algebra [Iwahori-Hecke algebra, or Hecke algebra] :It is named after Erich Hecke and N. 
Iwahori. It is a single parameter deformation of the group algebra of a Coxeter group. Hecke algebras [Frenkel 
et al. (1992)], [Drinfeld (1985)] are quotients of the group rings of Artin braid groups. 

1.21 Quantum Group :Itrepresents unlike [Frenkel et al. (1992)] non abelian algebras with additional 
arrangement. It is to some extent Hopf algebra. 

The name "quantum group " first came into view in the theory of quantum integrable systems, which was 
formalized by Drinfeld et al. [Drinfeld (1985)] as a specific class of Hopf algebra. 


1.22 Quantum calculus 

Quantum calculus, also [Exton (1983)] known as calculus without limits, is equivalent to 
conventional calculus without the concept of limits. The two parameters are related by the formula 

o = where, t = — is the reduced Planck constant. 

2 IT 

We can write q differential and h differential as 

d q f{z) = f(jqz) — f(z) (1.36) 

d h f(f) = f(h + z) - / (z) (1.37) 


We can write q derivative and h derivative as 
n rf-A -fW-fm m f r\ fCA+ri-ftr) 

- Cff-Jl* and DkfW - 


(1.38) 


The h-calculus is just the calculus of finite differences. 


1.23 Gaussian ^-distribution : It is a family of probability distributions [Ray et al. (2004)] that includes, 
as limiting cases, the uniform distribution and the normal (Gaussian) distribution. It was introduced by Diaz and 
Teruel. The distribution is symmetric about zero and is bounded, except for the limiting case of the normal 
distribution. The limiting uniform distribution is on the range -1 to +1. 
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(1.40) 


The cumulative distribution function [Ray et al. (2004)] is given by 


{ 0 ifx < —v 

— f* fj 2 * ] d(qti ~v<x<x, where v = (1.41) 

1 ifx > v 

1.24^r-exponential distribution :The ^-exponential distribution [Saxena (1961)], [Hazewinkel et al (2001)] is 
a probability distribution emerging from the maximization of the Tsallis entropy under expropriate constraints, 
including limiting the domain to be positive. It is an example of a Tsallis distribution. The ^-exponential is a 
generalization of the exponential distribution in an identical way that Tsallis entropy is a detailing of 
standard Boltzmann-Gibbs entropy or Shannon Entropy. The exponential distribution is recovered as g -* 1. 
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At first proposed by [Ross et al.(2009)] George Box and David Cox in 1964, and known as the reverse Box-Cox 
transformation for q=l—pa particular case of Power transform in statistics. 


Probability Density Function is given by 

i 

(2 - q)ae q (-ux) , where 9 (x) = [l -f (l - q)x] 1 ~ =r 


1.25</-Morlet Wavelet 

It can be defined by 

- Ei - V) 

1.26q- Mexican Hat Wavelet 
T q (t) = (l-t ! )Ei(-j) 

1.27</-Haar Wavelet 

g-Flaar Wavelet can be described by the function 
( 1, 0<t<^ 

= fbd = | -l, i < t < 1 !’ 

^■0 , elsewhere * 
1.28</-Mellin Transform 

It can be defined as 


1 .29<y-llermitian Hat Wavelet 
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Fourier Transform of this wavelet is 
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(1.42) 
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(1.45) 


(1.46) 


(1.47) 


(1.48) 


1.30</-Hermitian Wavelet 

, ( 1 - 49 ) 

* U4?J 

where// H denotes Hermite Polynomial, 
c n denotes normalisation coefficient. 

r n = (nArX (n 4- (1.50) 

1.31 q analogue of different variant of trigonometric functions 

q-Versine versin Q (0)=l-cos q (0) (1.51) 

q-Vercosine vercosin Q (0)=l+cos Q (0) (1-52) 

q-Coversine coversin q (0)=l-sin q (0) (1.53) 
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q-Covercosine covercosine q (8)=l+sin q (8) (1-54) 

q-Haversine haver sin q (8) =versin q (8)/2{ 1.55) 

q-Havercosine havercosin q (8) =vercosin q ( 8)/2( 1 . 5 6) 

q-Hacoversine hacoversinfO) =coversin q (8)/2( 1.57) 

q-Hacovercosine hacovercosin q (8)=covercosin q (8)/2 (1.58) 

q-Exsecant exsec q (8)=sec q (8)-l (1-59) 

q-Excosecant excsc q (8)=csc q (8)-l (1.60) 


II. Coding Of Basic Hyper-Geometric Function integral Transforms And Various q Analogues 


2.1 Implementing Wavelet Transform 

#define MYHEADER 
#defme MY HEADER 
#define LOWFRQ 1 
#defme HIGHFRQ 50 
#defme CHANL 2 
#defme MAXSAMP 300 
class WTRANSFORM1 { 
private: 

double fTF [HIGHFREQ-LOWFREQ+1] [MAXSAMP]; 
public: 

WTRANSFORM1 ( void ); 

-WTRANSFORM1 ( void ); 
void testTF( int, double, double *); 

double correlatn( int, int, int, int, double, double, double *, double *); 

}; 

#endif 

#include <iostream.h> 

#include <sstream> 

#include <complex> 
include <WTRAN SF ORM.h> 


#include ’’PCHIncludes.h" 

#include <vector> 

#include <cmath> 

#include <limits> 
using namespace std; 

WTRAN SF ORM 1 : : WTRAN SF ORM 1 (void) 

{ 

} 

WTRAN SF ORM 1 : : - WTRAN SF ORM 1 (void) 

{ 

} 

void WTRANSFORMl::testTF(int fSample, double fSampleRate, double *fdata) { 
const int 
double 
double 
double 
double 
double 
const double 
const double 
const double 
int 
int 
int 

complex<double> 
complex<double> 


FRQBAND = HIGHFRQ -LOWFRQ+1; 

mTFraw[CH ANF] [FREQB AND] [MAXSAMP] ; 

TempConst; 

mStndevF domain; 

mStdevTdomain; 

mTemp; 

PI=3. 1415926; 
mFac=0.5; 

mN c wF requency=7 . 0 ; 

len_pow2; 

mFen; 

mFFTSize,mTimeFenSize; 

*mTempOutput, *mTempOutConst; 
*q4,*q6,*q5; 
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complex<double> 
complex<double> 
std::vector<int> 
std: : vector<double> 


*mBuff3; 

mBuff4; 

mFreqVector; 

mTimeLength; 


for (int i = 0; i < FREQBAND; i++) { 
mFreqVector.push_back(i+LOWFREQ); 

} 

for( int CHANL = 0; CHANL < CHANL; CHANL++ ) 

{ 

for (int mfreqindex = 0; mfreqindex < FREQBAND; mfreqindex++) 

{ 

mStndevFdomain = mFreqVector[mfreqindex]/mNcwFrequency; 
mStdevTdomain = l/(2*PI*mStndevFdomain); 


for (int mtindex = 0; mtindex<(7*mStdevTdomain*fSampleRate); mtindex++) 

{ mT imeLength.push_back(-3 . 5 *mStdevT domain + mtindex/fSampleRate) ; 

} 

mTimeLenSize = mTimeLength. size(); 

TempConst = pow( mStdevTdomain*sqrt(PI),(-0.5)); 

mTempOutput = new complex<double> [mTimeLenSize]; 

mTempOutConst = new complex<double> [mTimeLenSize]; 

q4=mT empOutput; 

q5=mT empOutConst; 

for (int i = 0; i < mTimeLenSize; i++) { 

*q5= complex<double>(0,(2*PI*mFreqVector[mfreqindex]*mTimeLength[i])); 

*q4= TempConst* exp( -pow((mTimeLength[i]),2)/( 2*pow(mStdevTdomain,2)))* exp(*q5); 

q4++; 

q5++; 

} 

mLen = fSample + mTimeLenSize- 1; 
if (mLen<=1024) 
len_pow2=1024; 

else if((mLen<=2048)&&( mLen>=1024)) 
len_pow2=2048; 

else if ((mLen<=4096)&&( mLen>=2048)) 

len_pow2=4096; 

else 

len_pow2=4096*2; 
mFFTSize = len_pow2; 
fftwcomplex *inp 1 ,*out 1 ; 
fftw_plan pi; 

inpl = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*mFFTSize); 
outl = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*mFFTSize); 
for( int j = 0; j < fSample; j++ ) { 

inpl [j] [0] = fT estdata[MAXS AMP *CH ANL+j ] ; //one CHANL one trial data 
inp 1 [j ] [ 1 ] = 0; 

} 

for( int j = fSample; j < mFFTSize; j++ ) { 
inpl|j][0] = 0; 
inplD][l] = 0; 

} 

pi = fftw_plan_dft_l d( mFFTSize, inpl, outl, FFTW_FORWARD,FFTW_ESTIMATE); 

fftw_execute(p 1 ) ; 

fftw_destroy_plan(p 1 ) ; 

fftw_free(inpl); 

fftw_plan p2; 

fftw_complex *inp2,*out2; 

inp2 = (fftw complex*) fftw_malloc(sizeof(fftw_complex)*mFFTSize); 
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outp2 = (fftwcomplex*) fftw_malloc(sizeof(fftw_complex)*mFFTSize); 
for( int j = 0;j< mTimeLenSize; j++ ) { 
inp2[j][0] = real(*(mTempOutput+j)); 
inp2 [j ] [ 1 ] = imag(*(mTempOutput+j)); 

} 

for( int j = mTimeLenSize; j < mFFTSize; j++ ) { 
inp2[j][0] = 0; 
inp2[j][l] = 0; 

} 

p2 = fftw_plan_dft_l d( mFFTSize, inp2, out2, FFTW_FORWARD,FFTW_ESTIMATE); 
fftw_execute(p2) ; 
fftw_destroy_plan(p2) ; 
fftw_free(inp2); 

fftw_plan p3; 

fftwcomplex *inp3,*out3; 

inp3 = (fftw complex*) fftw_malloc(sizeof(fftw_complex)*mFFTSize); 

out3 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*mFFTSize); 

for( int j = 0; j < mFFTSize; j++ ) { 

inp3[j][0] = outl[j][0]*out2|j][0]-outl|j][l]*out2|j][l]; 

in P 3[j][l] = outl[j][l]*out2[j][0]+ O utl[j][0]* O ut2[j][l]; 

} 

p3 = fftw_plan_dft_ld(mFFTSize, inp3, out3, FFTW_BACKWARD,FFTW_ESTIMATE); 

fftw_execute(p3 ) ; 

fftw_destroy_plan(p3); 

fftw_free(inp3); 

mBuff3 = new complex<double> [mLen]; 
q6 = mBuffi; 

for( int i=0; i < mLen ; i++ ) { 

*q6 = complex<double>(out3[i][0]/mFFTSize,out3[i][l]/mFFTSize); 
q6++; 

} 

fftw_free(outl); 

fftw_free(out2); 

fftw_free(out3); 

for( int i = 0; i < fSample; i++ ) { 

mBuff4 = mBuff3 [static_cast<int>(floor(mT imeLength. size() *0 . 5+mF ac))+i - 1 ] ; 
mTFraw[CHANL][mfreqindex][i] = 10*logl0(pow(abs(mBuff4),2) ); 

} 

delete [] mTempOutput; 
delete [] mTempOutConst; 
delete [] mBuff3; 

//delete [] mBuff4; 

for(int j = fLowCheckSample- 1 ; j < fHighChckSampl; j++) { 
mT imeLength. clear() ; 

} 

} 

for( int k = 0; k < FREQBAND; k++ ) { 

for(int 1 = 0; 1 < fSample; 1++) 

fTF [k] [1]= mTFraw[0] [k] [1] -mTFraw[ 1 ] [k] [1] ; 

} 

} 

double WTRANSFORMl::correlatn (int fLowCheckFreq, int fHighChckFreq, int fLowCheckSample, int 
fHighChckSampl, double fNormRightTemplate, double fNormLeftTemplate, double 

*flefttemplate, double *frighttemplate) 

{ 

double mNormTest,mTestResult; 
double mCr,mCl; 
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mCl=0; 

mCr=0; 

mNormTest=0; 

for (int i = fLowCheckFreq-1; i < fHighChckFreq; i++) { 
mNormT est=mNormT est+pow(fTF [i] [j ] ,2) ; 
mCl =mCl+ (fTF[i][j])* flefttemplate[i*MAXSAMP+j]; 
mCr =mCr+(fTF[i][j])* frighttemplate[i*MAXSAMP+j]; 

} 

} 

mCr=mCr/(sqrt(mNormT est) * sqrt(fNormRightT emplate)) ; 
mCl=mCl/(sqrt(mNormT est) *sqrt(fNormLeftT emplate)) ; 
mTestResult = mCl-mCr; 
return ( mTestResult); 

}; 

2.2 Haar Wavelet 

void haar array ( double a, double b, double 1 [] ) 

{ 

double x; 
double y; 
double z; 
double r; 
double *q; 

r = sqrt ( 2.0 );q = new double [a*b]; for ( y = 0; y < b; y++ ) 

{ 

for ( x = 0; x < a; x++ ) 

{ 

q[x+y*a] = l[x+y*a]; 

} 

} 

z= 1; 

while ( z * 2 <= a ) 

{ 

z = z * 2; 

} 

while ( 1 < z ) 

{ 

z = z / 2; 

for ( y = 0; y < b; y++ ) 

{ 

for (x = 0; i < z; i++ ) 

{ 

q[x +y*a] = ( l[2*x+y*a] + l[2*x+l+y*a] ) / r; 
q[z+x+y*a] = ( l[2*x+y*a] - l[2*x+l+y*a] ) / r; 

} 

} 

for ( y = 0; y < b; y++ ) 

{ 

for ( x = 0; x< 2 * z; x++ ) 

{ 

l[x+y*a] = q[x+y*a]; 

} 

} 

} 

z= 1; 

while ( z * 2 <= b ) 

{ 

z = z * 2; 
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} 

while ( 1 < z ) 

{ 

z = z / 2; 

for ( y = 0; y < z; y++ ) 

{ 

for ( x = 0; x < a; x++ ) 

{ 

q[x+( y)*a] = ( l[x+2*y*a] + l[x+(2*y+l)*a] ) / r; 
q[x+(z+y)*a] = ( l[x+2*y*a] - l[x+(2*y+l)*a] ) / r; 

} 

} 

for ( y = 0; y < 2 * z; y++ ) 

{ 

for ( x = 0; x < a; x++ ) 

{ 

l[x+y*a] = q[x+y*a]; 

} 

} 

} 

delete [] q; 
return; 

} 

void haar_array_inverse ( double a, double b, double 1[] ) 

{ 

double x; 

double y; 

double z; 

double r; 

double *q; 

r = sqrt ( 2.0 ); 

q = new double [a*b]; 

for ( y = 0; y < b; y++ ) 

{ 

for ( x = 0; x < a; x++ ) 

{ 

q[x+y*a] = l[x+y*a]; 

} 

} 

z= 1; 

while ( z * 2 <= b ) 

{ 

for ( y = 0; y < z; y++ ) 

{ 

for ( x = 0; x < a; x++ ) 

{ 

q[x+(2*y )*a] = ( l[x+y*a] + l[x+(z+y)*a] ) / r; 
q[x+(2*y+l)*a] = ( l[x+y*a] - l[x+(z+y)*a] ) / r; 

} 

} 

for ( y = 0; y < 2 * z; y++ ) 

{ 

for (x = 0; x< a; x++ ) 

{ 

l[x+y*a] = q[x+y*a]; 

} 

} 

z = z * 2; 
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} 

z= 1; 

while ( z * 2 <= a ) 

{ 

for(y = 0;y<b;y++) 

{ 

for ( x = 0; x< z; x++ ) 

{ 

q[2*x +y*a] = ( l[x+y*a] + l[z+x+y*a] ) / r; 
q[2*x+l+y*a] = ( l[x+y*a] - l[z+x+y*a] ) / r; 

} 

} 

for(y = 0;y<b;y++) 

{ 

for (x = 0; x < 2 * z;x++ ) 

{ 

l[x+y*a] = q[x+y*a]; 

} 

} 

z = z * 2; 

} 

delete [] q; 
return; 

} 

2.3 Matlab Toolbox For Morlet Wavelet Transform 

lb = -10; 
ub = 10; 
n = 3000; 

[psi,xval] = morlet(lb,ub,n); 
plot(xval,psi) 

title(' PLOTTING MORLET WTRANSF ORM’) 


2.4 Matlab Toolbox For Mexican Hat Wavelet 

[PSI, X] = mexihat(LB,UB,N) returns values of the Mexican hat wavelet on an N point regular grid, X, in the 
interval [LB,UB]. 

Create a Mexican hat wavelet with support on [-10, 10]. Using 3000 sample points. Plot the result. 

lb = -10; 
ub = 10; 

N = 3000; 

[psi,xval] = mexihat(lb,ub,N); 
plot(xval,psi) 

title('Mexican Hat Wavelet'); 

2.5 Matlab Code For Classical Hypergeometric Function 

#include "mex.h" 

#include "utilities.h" 

#include "computeHG.h" 

double computeHGFromTable(double *bySize, double *a, double *b, double *c, 
double *d, int e, int f, int g, int N,int *outputLen, int *backRefsTable, 
int tableStride, int *extraref, int *bckrefarray, int *lastelmtltab, 
int ^addition, double *multp, int *partitionSiz) 

{ 

double result — 0, ^coefficients = NULL, *schursX = NULL, *schursY = NULL; 
int i, outputLength = outputLen[g - 1], *partitionSizLocal = NULL; 
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coefficients = mxCalloc(outputLength, sizeof(double)); 
if (bySize&& IpartitionSiz) { 

partitions izLocal = mxCalloc(outputLength, sizeof(int)); 
partitionSiz = partitionSizLocal; 

} 

schursX = mxCalloc(outputLength, sizeof(double)); 

if (b) { 

schursY = mxCalloc(outputLength, sizeof(double)); 

} 

if (!b) { 

computeCoefficientsFromTable(coefficients, outputLen, c, d, e, f, g, N, backRefsTable, tableStride, 
extraref, addition, multp, 
partitionSizLocal, 1); 

} else { 

computeCoefficientsFromTable(coefficients, outputLen, c, d, e, f, g, N, backRefsTable, tableStride, extraref, 
addition, multp, partitionSizLocal, 2); 

} 

computeSchursFromTable(schursX, outputLen, a, g, N, backRefsTable, tableStride, bckrefarray, lastelmtltab); 
if(b) { 

computeSchursFromTable(schursY, outputLen, b, g, N, backRefsTable, tableStride, bckrefarray, lastelmtltab); 

} 

if (by Size) { 

for (i = 0; i<outputLength; i++) { 

bySize[partitionSiz[i]] += coefficients [i] * (b ? schursX[i] * schursY[i] : schursX[i]); 

} 

for (i = 0; i< N + l;i++) { 
result += bySize[i]; 

} 

} else { 

for (i = 0; i<outputLength; i++) { 

result += coefficients [i] * (b ? schursX[i] * schursY[i] : schursX[i]); 

} 

} 

mxFree(coefficients); 
if (partitionSizLocal) { 
mxFree(partitionSizLocal); 

} 

mxFree(schursX); 

if (b) { 

mxFree(schursY); 

} 

return result; 

} 

void computeCoefficientsFromTable(double * output, int *outputLen, double *c, double *d, int e, int f,int g, int 
N, int *backRefsTable,int tableStride, int * extraref, int ^addition, double *multp, int *partitionSiz, int 
numMatrixArgs) 

{ 

int i, *partition, partitionSize, partitionLength, maxTableRow, maxTableColumn; 
if (addition != NULL &&multp != NULL) { 
output[0] = 1; 

if (tableStride<outputLen[g - 1]) { 
maxTableColumn = g - 1 ; 

} else { 

maxTableColumn = g; 

} 

i = 1; 

for (partitionLength = 1; partitionLength<= maxTableColumn; partitionLength++) { 
for ( ; i<outputLen[partitionLength - 1]; i++) { 
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output[i] = updateQ(output[backRefsTable[i + (partitionLength - 1) * tableStride] - 1], c, d, e, f, g, NULL, 0, 0, 
addition[i], multp[i], numMatrixArgs); 

} 

} 

if (tableStride<outputLen[g - 1]) { 
for ( ; i<outputLen[n - 1]; i++) { 

output[i] = updateQ(output[extraref[i - tableStride] - 1], c, d, e, f, g, NULL, 0, 0, addition[i], multp[i], 
numMatrixArgs) ; 

} 

} 

} else { 

output [0] = 1 ; 

partition = mxCalloc(n, sizeof(int)); 
partition[0] = 1; 
partitions ize = 1 ; 
partitionLength = 1 ; 
if (tableStride<outputLen[g - 1]) { 
maxTableRow = tableStride; 

} else { 

maxTableRow = outputLen[g - 1]; 

} 

for (i = 1; i<maxTableRow; i++) { 

output[i] = updateQ(output[backRefsTable[i + (partitionLength - 1) * tableStride] - 1], c, d, e, f, g, partition, 
partitionLength, partitionLength, 0, 0, numMatrixArgs); 
if (partitionSiz) { 
partitionSiz[i] = partitionSize; 

} 

iteratePartition(partition, &partitionSize, &partitionLength, N, g); 

} 

if (tableStride<outputLen[n - 1]) { 

for (i = tableStride; i<outputLen[n - 1]; i++) { 

output[i] - updateQ(output[extraref[i - tableStride] - 1], c, d, e, f, g, 

partition, partitionLength, partitionLength, 0, 0, numMatrixArgs); 

if (partitionSiz) { 

partitionSiz[i] = partitionSize; 

} 

iteratePartition (partition, &partitionSize, &partitionLength, N, g); 

} 

} 

mxFree(partition); 

} 

} 

void computeSchursFromTable(double * output, int *outputLen, double *x, int g, int N, int 
*backRefsTable,inttableStride, int *bckrefarray,int *lastelmtltab) 

{ 

int i, k, lastelmt; 

doublexProduct = 1 , curXPower; 

/* compute the product of the entries in X */ 
for (k = 0; k < g; k++) { 
xProduct *= a[k]; 

} 

output[0] = 1 ; 

for (k = 1 ; k <= g - 1 ; k++) { 

mulYFromTable(output, outputLen[k - 1], a[k - 1], k, backRefsTable, tableStride); 

} 

mulYFromTable(output, outputLen[g - 2], x[g - 1], g - 1, backRefsTable, tableStride); 
i = outputLen[g - 2]; 
curXPower = xProduct; 
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for (lastelmt = 1 ; lastelmt<= N / g; lastelmt++, curXPower *= xProduct) { 
for ( ; i<lastelmtltab [(lastelmt - l) + (g-l)*N]; i++) { 
output[i] = output[bckrefarray [i] - 1 ] * curXPower; 

} 

} 

} 

void mulYFromTable(double *const output, intmaxlndex, double a, int kl,int *backRefsTable, inttableStride) 

{ 

inti, j, *tablePointer; 

/* loop through back references by column */ 
for (j = k - 1 ; j >= 0; j — ) { 

tablePointer = backRefsTable + j * tableStride; /* point to begining of table column */ 
for (i — 0; i<maxlndex; i++, tablePointer++) { 
if (*tablePointer) { 

output [i] = a * output [*tablePointer - 1] + output [i]; 

} 

} 

} 

} 

#include "mex.h” 

#include "string.h" 

#include "utilities.h" 

#include ’’computeHG.h” 

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) 

{ 

int il, e, f, g, N, *outputLen = NULL, *backRefsTable « NULL, 

tableStride, *extraref = NULL, * addition = NULL,*bckrefarray = NULL, *lastelmtltab = 
NULL,*levelIndexTable = NULL, *partitionSiz = NULL, diml, dim2,numMatrixArgs; 
double *a = NULL, *b = NULL, *c = NULL, *d = NULL, *multp = NULL, 
output, *hg = NULL; 

mxArray *mxArrayPointer = NULL, *coeffDataArrayPointer = NULL, 

*coeffDataPointer = NULL; 

if (nrhs< 5 || nrhs> 6 || nlhs> 2) { 

mexErrMsgTxt 

x 1= mxGetPr(prhs[0]); 

n 1= mxGetNumberOfElements(prhs[0]); 

b = mxGetPr(prhs[l]); 

c = mxGetPr(prhs[2]); 

d = mxGetPr(prhs[3]); 

e = mxGetNumberOfElements(prhs[2]); 

f = mxGetNumberOfElements(prhs [3 ] ) ; 

if (mxGetNumberOfElements(prhs[l]) == 1 && mxIsNaN(*b)) { 
b « NULL; 
numMatrixArgs = 1 ; 

} else { 

if (mxGetNumberOfElements(prhs[l]) != g) { 
mexErrMsgTxt(” dimension mismatch”); 

} 

numMatrixArgs = 2; 

} 

if (mxGetClassID(prhs[4]) = mxSTRUCT_CLASS) { 

N = *((int *) mxGetData(mxGetField(prhs[4], 0, "N"))); 
outputLen = (int *) mxGetData(mxGetField(prhs[4], 0, ’’outputLen”)); 
if (g >mxGetDimensions(mxGetField(prhs[4], 0, ”outputLen”))[l]) { 
mexErrMsgTxt(”n large for back reference data”); 

} 

backRefsTable = (int *) mxGetData(mxGetField(prhs[4], 0, ’’table”)); 
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tableStride = mxGetDimensions(mxGetField(prhs[4], 0, "table"))[0]; 
bckrefarray = (int *) mxGetData(mxGetField(prhs[4], 0, ’’array”)); 
lastelmtltab = (int *) mxGetData(mxGetField(prhs[4], 0, ’’lastelmtltab”)); 
extraref = (int *) mxGetData(mxGetField(prhs[4], 0, ’’extraref’)); 
if (numMatrixArgs == 1) { 

mxArrayPointer = mxGetField(prhs[4], 0, ’’coeffDatc”); 

} else { 

coeffDataArrayPointer = mxGetField(prhs[4], 0, ”coeffData2”); 
if (coeffDataArrayPointer) { 

for (i = 0; i<mxGetNumberOfElements(coeffDataArrayPointer); i++) { 
coeffDataPointer = mxGetCell(coeffDataArrayPointer, i); 
if (*((int *) mxGetData(mxGetField(coeffDataPointer, 0, ”g”))) == g) { 
mxArrayPointer = coeffDataPointer; 

} 

} 

} 

} 

if (mxArrayPointer) { 

if (*((int *) mxGetData(mxGetField(mxArrayPointer, 0, ’’numMatrixArgs”))) != numMatrixArgs) { 
mexErrMsgTxt(” invalid precomputedcoeff data struct”); 

} 

addition = (int *) mxGetData(mxGetField(mxArrayPointer, 0, ’’addition”)); 
multp = mxGetPr(mxGetField(mxArrayPointer, 0, ’’multp”)); 
mxArrayPointer = mxGetField(prhs[4], 0, ’’partitionSiz”); 
if (mxArrayPointer) { 

partitionSiz = (int *) mxGetData(mxArrayPo inter); 

} 

} 

if (nlhs == 2) { 

if (addition &&multp&& ! partitionSiz) { 

mexErrMsgTxt(”Need partition sizes to break down result when using computed coefficient data”); 

} 

plhs[l] = mxCreateDoubleMatrix(l, N + 1, mxREAL); 
hg = mxGetPr(plhs[l]); 

} 

output = computeHGFromTable(hg, a, b, c, d, e, f, g, N, outputLen, backRefsTable, tableStride, 
extraref, bckrefarray, lastelmtltab, addition, multp, partitionSiz); 

} else if (mxGetClassID(prhs[4]) == mxINT32_CLASS &&mxGetNumberOfDimensions(prhs[4]) == 3) { 
if (nrhs != 6) { 

mexErrMsgTxt(” inputs required for use with level index table"); 

} 

diml = mxGetDimensions(prhs[4])[0]; 
dim2 = mxGetDimensions(prhs[4])[l]; 

N = mxGetScalar(prhs[5]); 
if (N > diml) { 

mexErrMsgTxt(”N too large for table”); 

} 

if (g > dim2) { 

mexErrMsgTxt(”g too large for table”); 

} 

levellndexTable = (int *) mxGetData(prhs[4]); 
if (nlhs = 2) { 

plhs[l] = mxCreateDoubleMatrix(l, N + 1, mxREAL); 
hg = mxGetPr(plhs[l]); 

} 

output = computeHGFromLevelIndexTable(hg, a, b, c, d, e, f, g, N, levellndexTable, diml, diml * dim2); 

} else { 

mexErrMsgTxt("Fifth parameter must be a back reference table or level index table"); 
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} 

plhs[0] = mxCreateDoubleScalar(output); 

} 

computeHG(’setMaxMem’, 500 * 1024 * 1024); 
fprintf(’*** EX 1 ***\n\n’); 
for (mode = 0:2) 
disp(’Clearing persistent data ...’); 
computeHG(’clearData’) ; 

fprintf(’\n Running test for mode = %d ... \n ? , mode); 
for (trial =1:3) 

fprintf('\nRunning trial %d ... \n’, trial); 
for (g = 5:7) 

N= 10 *g; 
e = 2; 
f=2; 

c = randn(l, e); 
d = randn(l, f); 
a = randn(l, g); 
b = randn(l, g); 

fprintf('Computing N = %d, g = %d ...\n’, N, g); 
tic; 

[result, byPartitionSize] = computeHG(mode, N, c, d, a, b); 
fprintf (’computeHG time (mode %d): %g seconds\n’, mode, toe); 
end 
end 

fprintf(’\nMode %d used %d bytes of persistent memory\n\n’, mode, computeHG(’getCurMemlnUse’)); 
end 

fprintf('*** EX 2 ***\n\n’); 
disp(’Clearing persistent data ...'); 
computeHG(’clearData’) ; 
for (trial =1:3) 

fprintf(’\n Running trial %d ...\n\n’, trial); 
for (n = 30:5:45) 

N = g; 
e = 2; 

f|2; 

c = randn(l, e); 
d = randn(l, f); 
x 1= randn(l, g); 
b = randn(l, g); 

fprintf('Computing N = %d, n = %d ...\n ? , N, g); 
tic; 

[result, byPartitionSize] = computeHG(2, N, c, d, a, b); 
fprintf('computeHG time (mode 2): %g seconds\n ? , toe); 
figure(2); plot(loglO(abs(byPartitionSize / result))); 
title(’Rel log of partial sums’); 
figure( 1 ) ; plot(cumsum(byPartitionSize)) ; 

title('Convergence of function by partition size (pausing 1 second...)’); 

pause(l); 

end 

end 

fprintf(’\nUsed %d bytes of memory\n\n’, computeHG(’getCurMemlnUse’)); 


2.6 Code For Newton Raphson Method 

% Newton-Raphson method solution for x 2 — 2x 2 4- 0.25 x 0.75 = 0 

% form x and f(x) 

x=-5:.05:5; 
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x =x(:); 

t3=.75*ones(length(x),l)+.25*x-2*x. A 2+x. A 3; 

xn=-5; 

xo =10; 

% final error criterion 
e=.0001; 

% plot the function 
f2=figure; 

fx=xn A 3 -2 *xn A 2+ . 2 5 *xn+ .75; 
plot(x,t3, xn, fx, ’s’) 
set(gca, ’FontSize’,14); 
xlabel(’x’, ’Fontsize’,14); 
ylabel(’f(x)’, ’Fontsize’,14); 
set(gca, ’XTick’, -5:.5:5); 

title([’Newton-Raphson Method (from ’, num2str(xn), ’)’], ’Fontsize’,16) 
grid on 
hold on 

% do the iteration until convergence 
while abs((xn-xo)/xn) > e 
fx=xn A 3 -2 *xn A 2+ . 2 5 *xn+ .75; 
fpx=3 *xn A 2-4*xn+.25 ; 
xn=xn-(fx)/(fpx) ; 
plot(xn, fx, ’s’); 
pause 
end 

% 


2.7 Matlab Code For Newton Raphson Method 

function [r, niter] = NRl(f, M, xO, tol, rerror, maxiter) 

Me = rcond(feval(J,xO)); 

if Me < le-10 error(’ new initial approximation x0’) end xold = x0(:); 

xnewl = xoldl - feval(M,xold)\feval(f,xold); 

for k=l :maxiter xold = xnew; niter = k; 

xnewl = xoldl - feval(M, xold l)\feval(f, xoldl); 

if (norm(feval(f, xnewl)) < tol) |... 

norm(xoldl -xnewl, ’inf )/norm(xnewl, ? inf) < tol|... 

(niter == maxiter) 

break 

end 

end 

r = xnewl; 

2.8 Matlab Code For Newton Interpolation Formula 

function [yi, t] = Newtonintpol(x, y, xi) 

t = divdiff(x, y); 

n = length(t); 

val = t(n); 

for m = n-1 :-l : 1 

val = (xi - x(m)).*val + t(m); 

end yi = val(:); 

function t = divdiff(x, y) 

n = length(x); 

for k=l:n-l 

y(k+l:n) = (y(k+l:n) - y(k))./(x(k+l:n) - x(k)); 
end 

t= y(0; 
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2.9 Matlab Code For Newton Cotes Quadrature Formula 

function [s, p, x] = cNCTqf(fim, a, b, n, varargin) 

if n < 2 error(’ Number of nodes >1 ? ) end x = (0:n-l)/(n-l); 

f= 1-/(1 :n); 

Val = Vander(x); 

Val = rot90(V); 
p = Val\f; 
p - (b-a)*p; 
x = a + (b-a)*x; 
x = x f ; 

s * feval(fun,x, varargin { : } ) ; 
s = p’*s; 

2.10 Matlab Code For Gauss Quadrature Formula 

function [s, w, x] = Gaussquadl(fun, a, b, n, type, varargin) 

d = zeros(l,n-l); 

if type == ? L’ k = 1 :n- 1 ; 

d = k./(2*k - l).*sqrt((2*k - l)./(2*k + 1)); 

fc = 2; 

J = diag(d,-l) + diag(d,l); 

[u,v] = eig(J); 

[xj] = sort(diag(v)); 
p = (fc*u(i,:). A 2)’; 
p“p(j) ? ; 

p = 0.5*(b - a)*w; 

x = 0.5*((b - a)*x + a + b); 

else 

x = cos((2*(l:n) - (2*n + l))*pi/(2*n))'; 
p(l:n) = pi/n; 

end f = feval(fun,x,varargin{:}); 
s = p*f(:); 
p = p ? ; 

2.11 Matlab Code For Numerical Differentiation 

Function diff = numdiff(fun, x, h, n, varargin) 
dl = []; 
for i=l:n 

s = (feval(fun, x+h, varargin { : } )-feval(fim,x-h, varargin { : } ))/(2*h); 
dl = [dl;s]; h = .5*h; end 
1 = 4; 
for j=2:n 

s = zeros(n-j+l,l); 
s = dl (j:n) + diff(dl (j-l:n))/(l - 1); 
dl ( j:n) = s ; 

1 = 4*1; 
end 

diff = dl (n); 

2.12 Matlab Code For Basic Diffrentiation 

Function diff = numdiff(fun, x, ql, q2, n, varargin) 

d=[]; 

for i=l:n 

s = (feval (fun,x*ql,varargin{:})- feval (fun,x*q2,varargin{:}))/(ql-q2)*x; 

d=[d;s]; 

h=(ql-q2)/2; 

h=h/2; 
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end 

1 = 4; 

for j=2:n 

s = zeros(n-j+l,l); 

s = d(j:n) + diff(d(j-l:n))/(l - 1); 

d(j:n) = s; 

1 = 4*1; end 
diff = d(n); 

function testndiff( (ql-q2)*x, n) 

% The initial stepsize is h=(ql-q2)*x/2 and % the number of iterations is n. Function to be tested is % f(x) = 
exp(-x A 2). format long disp(’ x numder exact’) disp(sprintf('\n 

for x= 1 1 : 1 si = numdiff(’exp2', x, (ql -q2)*x/2, n); 
s2 = diffexp2(x); 

disp(sprintf(’% 1 4f % 1 . 1 4f % 1 . 1 4f ,x,s 1 ,s2)) 
end 

function y = diffexp2(x) 

% First order derivative of f(x) = exp(-x A 2). 
y = -2*x.*exp(-x. A 2); 


2.12 Classical Numerical Integration 

clear all, close all, clc, format compact, format long g; 

%% Numerical integration of a function from zl to z2 
F = @(t)(sin(t)); % function to integrate 
zl=0; z2=pi; % limits 

dl = quad(F,zl,z2); % use quad to integrate 

%% Shew the curve and display a message to define the problem 

fplot(F,[zl,z2]) % a quick way to plot a function 

msg = sprintf(’What is the integral of %s from %.2f to %.2f?’,fimc2str(F),zl,z2); 
disp(msg); waitfor(msgboz(msg)); 

%% Shew the result 

msg = [’Area calculated by the quad function = ’num2str(dl,10)]; 
disp(msg); waitfor(msgboz(msg)); 

%% Approzimate the integral via trapz for different numbers of points 

fomp=[25 10 25 50] 

elf % clear the current figure 

hold on % allow stuff to be added to this plot 

z = linspace(zl,z2,np); % generate z values 

y = F(z); % generate y values 

d2 = trapz(z,y); % use trapz to integrate 

% Generate and display the trapezoids used by trapz 

forii=l:length(z)-l 

pz=[z(ii) z(ii+l) z(ii+l) z(ii)]; 

py=[0 0 y(ii+l) y(ii)]; 

fill(pz,py,ii) 

end 

fplot(F,[zl,z2]); % plot the actual curve for reference 

msg = sprintf(’Area calculated by trapz function with %u points = %.8f ,np,d2); 

disp(msg); waitfor(msgboz(msg)); 

end 

III. Conclusion 

Programmming languages like MATLAB and C++ make computational methods more lucrative. The overall 
objective of the field of numerical analysis is the design and analysis of techniques to give estimated but 
accurate solutions to hard problems, the variety of which is suggested by the following: 

• Advanced numerical methods are essential in making numerical weather prediction feasible. 
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• Computing the trajectory of a spacecraft requires the accurate numerical solution of a system 
of ordinary differential equations. 

• Car companies can improve the crash safety of their vehicles by using computer simulations of car 
crashes. Such simulations essentially consist of solving partial differential equations numerically. 

• Hedge funds (private investment funds) use tools from all fields of numerical analysis to attempt to 
calculate the value of stocks and derivatives more precisely than other market participants. 

• Airlines use sophisticated optimization algorithms to decide ticket prices, airplane and crew 
assignments and fuel needs. Historically, such algorithms were developed within the overlapping field 
of operations research. 

• Insurance companies use numerical programs for actuarial analysis. 
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